﻿#pragma once

#include <learn/geom/geom_basis.h>

namespace xi
{

namespace geom
{

template <std::size_t N, bool Rat = false> class Geom_BezierSurface : public Geom_BasisSurface<N>
{
    static constexpr std::size_t D = N + Rat;
    using VecR = VecT<D, double>; // (wP,w)
    using Vec = VecT<N, double>;  // P

  public:
    static constexpr bool Rational = Rat;
    Geom_BezierSurface(Geom_BezierSurface &&) = delete;

    /**
     * @brief 基本构造函数
     *
     * @param[in] cpts
     * @param[in] uDeg
     * @param[in] vDeg
     * @param[in] normalize
     */
    Geom_BezierSurface(std::vector<VecR> &&cpts, std::size_t uDeg, std::size_t vDeg, bool normalize = true)
        : m_uDeg(uDeg), m_vDeg(vDeg), m_cpts(std::move(cpts))
    {
        ASSERT_INFO(uDeg == u_degree() && vDeg == v_degree(), "Degree not match.");
        ASSERT_INFO(m_cpts.size() == u_pole_num() * v_pole_num(), "Control points not match.");

        // 如果为有理 Bezier 曲面，则对点的坐标进行规范化 (P,w) -> (wP,w)
        if constexpr (Rat)
        {
            if (normalize)
            {
                for (auto &p : m_cpts)
                {
                    double w = p[N];
                    p[N] = 1.0;
                    p *= w;
                }
            }
        }
    }

    /**
     * @brief u 方向的次数
     *
     * @return std::size_t
     */
    std::size_t u_degree() const
    {
        return m_uDeg;
    }

    /**
     * @brief v 方向的次数
     *
     * @return std::size_t
     */
    std::size_t v_degree() const
    {
        return m_vDeg;
    }

    /**
     * @brief u 控制点个数
     *
     * @return std::size_t
     */
    std::size_t u_pole_num() const override
    {
        return m_uDeg + 1;
    }

    /**
     * @brief v 控制点个数
     *
     * @return std::size_t
     */
    std::size_t v_pole_num() const override
    {
        return m_vDeg + 1;
    }

    /**
     * @brief 获得控制点
     *
     * @param[in] u
     * @param[in] v
     * @return Vec
     */
    Vec pole(std::size_t u, std::size_t v) const override
    {
        if constexpr (Rat)
            return homogenous_to_cartesian(m_cpts[u * v_pole_num() + v]);
        else
            return m_cpts[u * v_pole_num() + v];
    }

    /**
     * @brief 设置控制点
     *
     * @param[in] u
     * @param[in] v
     * @param[in] p
     */
    void set_pole(std::size_t u, std::size_t v, const Vec &p) override
    {
        if constexpr (Rat)
            m_cpts[u * v_pole_num() + v] = cartesian_to_homogenous(p) * weight(u, v);
        else
            m_cpts[u * v_pole_num() + v] = p;
    }

    /**
     * @brief 获得权重
     *
     * @param[in] u
     * @param[in] v
     * @return double
     */
    double weight(std::size_t u, std::size_t v) const
    {
        if constexpr (Rat)
            return m_cpts[u * v_pole_num() + v][N];
        else
            return 1.0;
    }

    /**
     * @brief 设置权重
     *
     * @param[in] u
     * @param[in] v
     * @param[in] w
     */
    void set_weight(std::size_t u, std::size_t v, double w)
    {
        static_assert(Rat, "Rational surface not supported.");
        m_cpts[u * v_pole_num() + v] *= w / weight(u, v);
    }

    /**
     * @brief 获得无穷远点
     *
     * @param[in] u
     * @param[in] v
     * @return VecR
     */
    VecR inf_pole(std::size_t u, std::size_t v) const
    {
        static_assert(Rat, "Rational surface required.");
        return m_cpts[u * v_pole_num() + v];
    }

    /**
     * @brief 设置无穷远点
     *
     * @param[in] u
     * @param[in] v
     * @param[in] p
     */
    void set_inf_pole(std::size_t u, std::size_t v, const VecR &p)
    {
        static_assert(Rat, "Rational surface required.");
        m_cpts[u * v_pole_num() + v] = p;
    }

    /**
     * @brief 获得控制点
     *
     * @return const std::vector<VecR>&
     */
    const std::vector<VecR> &poles() const
    {
        return m_cpts;
    }

    /**
     * @brief 移动控制点
     *
     * @return std::vector<VecR>
     */
    std::vector<VecR> move_poles()
    {
        return std::move(m_cpts);
    }

    /**
     * @brief 参数范围
     *
     * @return std::array<double, 4>
     */
    std::array<double, 4> bounds() const override
    {
        return {0.0, 1.0, 0.0, 1.0};
    }

    /**
     * @brief 计算曲面上点的位置
     *
     * @param[in] u
     * @param[in] v
     * @return Vec
     */
    Vec d0(double u, double v) const override
    {
        auto p = de_casteljau2(m_cpts.data(), u_pole_num(), v_pole_num(), u, v);
        if constexpr (Rat)
            return homogenous_to_cartesian(p);
        else
            return p;
    }

    /**
     * @brief 计算曲面上点的切向量
     *
     * @param[in] u
     * @param[in] v
     * @return std::array<Vec, 3>
     */
    std::array<Vec, 3> d1(double u, double v) const override
    {
        std::size_t ud = u_degree();
        std::size_t vd = v_degree();
        std::vector<double> uKnots((ud + 1) * 2);
        std::vector<double> vKnots((vd + 1) * 2);
        std::fill(uKnots.begin(), uKnots.begin() + ud + 1, 0.0);
        std::fill(uKnots.begin() + ud + 1, uKnots.end(), 1.0);
        std::fill(vKnots.begin(), vKnots.begin() + vd + 1, 0.0);
        std::fill(vKnots.begin() + vd + 1, vKnots.end(), 1.0);
        if constexpr (Rat)
        {
            auto ders = rat_surface_deriv(m_cpts.data(), uKnots.data(), vKnots.data(), uKnots.size(), vKnots.size(),
                                          u_degree(), v_degree(), 1, u, v);
            return {ders[0][0], ders[1][0], ders[0][1]};
        }
        else
        {
            auto ders = surface_deriv_alg1(m_cpts.data(), uKnots.data(), vKnots.data(), uKnots.size(), vKnots.size(),
                                           u_degree(), v_degree(), 1, u, v);
            return {ders[0][0], ders[1][0], ders[0][1]};
        }
    }

    /**
     * @brief 计算曲面上点的二阶导数
     *
     * @param[in] u
     * @param[in] v
     * @return std::array<Vec, 6>
     */
    std::array<Vec, 6> d2(double u, double v) const override
    {
        std::size_t ud = u_degree();
        std::size_t vd = v_degree();
        std::vector<double> uKnots((ud + 1) * 2);
        std::vector<double> vKnots((vd + 1) * 2);
        std::fill(uKnots.begin(), uKnots.begin() + ud + 1, 0.0);
        std::fill(uKnots.begin() + ud + 1, uKnots.end(), 1.0);
        std::fill(vKnots.begin(), vKnots.begin() + vd + 1, 0.0);
        std::fill(vKnots.begin() + vd + 1, vKnots.end(), 1.0);
        if constexpr (Rat)
        {
            auto ders = rat_surface_deriv(m_cpts.data(), uKnots.data(), vKnots.data(), uKnots.size(), vKnots.size(),
                                          u_degree(), v_degree(), 2, u, v);
            return {ders[0][0], ders[1][0], ders[0][1], ders[2][0], ders[0][2], ders[1][1]};
        }
        else
        {
            auto ders = surface_deriv_alg1(m_cpts.data(), uKnots.data(), vKnots.data(), uKnots.size(), vKnots.size(),
                                           u_degree(), v_degree(), 2, u, v);
            return {ders[0][0], ders[1][0], ders[0][1], ders[2][0], ders[0][2], ders[1][1]};
        }
    }

    /**
     * @brief 序列化
     *
     * @return utils::Json
     */
    utils::Json dump_json() const override
    {
        utils::Json json;
        json["type"] = "bezier_surface";
        json["rational"] = Rat;
        json["uDegree"] = int(u_degree());
        json["vDegree"] = int(v_degree());
        json["poles"] = m_cpts;
        return json;
    }

  protected:
    std::size_t m_uDeg;
    std::size_t m_vDeg;
    std::vector<VecR> m_cpts;
};

} // namespace geom

} // namespace xi
